# Loading the csv files from dropbox
hamilton_final_price <- read.table("~/Dropbox/Replication Workshop/Student Papers/Saleuddin/Data/OriginalData/hamilton_final_price.csv", quote="\"")
hamilton_final_CornS <- read.table("~/Dropbox/Replication Workshop/Student Papers/Saleuddin/Data/OriginalData/hamilton_final_CornS.csv", quote="\"")
hamilton_final_CornF <- read.table("~/Dropbox/Replication Workshop/Student Papers/Saleuddin/Data/OriginalData/hamilton_final_CornF.csv", quote="\"")
hamilton_final_OatsS <- read.table("~/Dropbox/Replication Workshop/Student Papers/Saleuddin/Data/OriginalData/hamilton_final_OatsS.csv", quote="\"")
hamilton_final_OatsF <- read.table("~/Dropbox/Replication Workshop/Student Papers/Saleuddin/Data/OriginalData/hamilton_final_OatsF.csv", quote="\"")
hamilton_final_RyeS <- read.table("~/Dropbox/Replication Workshop/Student Papers/Saleuddin/Data/OriginalData/hamilton_final_RyeS.csv", quote="\"")
hamilton_final_LardS <- read.table("~/Dropbox/Replication Workshop/Student Papers/Saleuddin/Data/OriginalData/hamilton_final_LardS.csv", quote="\"")
hamilton_final_LardF <- read.table("~/Dropbox/Replication Workshop/Student Papers/Saleuddin/Data/OriginalData/hamilton_final_LardF.csv", quote="\"")
hamilton_final_TBill <- read.table("~/Dropbox/Replication Workshop/Student Papers/Saleuddin/Data/OriginalData/hamilton_final_TBill.csv", quote="\"")
hamilton_final_RyeF <- read.table("~/Dropbox/Replication Workshop/Student Papers/Saleuddin/Data/OriginalData/hamilton_final_RyeF.csv", quote="\"")
hamilton_final_RyeS <- hamilton_final_RyeS[-58]
hamilton_final_RyeF <- hamilton_final_RyeF[-58]
# converting them into time series objects
Price.ts = ts(data=hamilton_final_price, frequency = 3, start=c(1920,2), end=c(1939,2)) 
CornS.ts = ts(data=hamilton_final_CornS, frequency = 3, start=c(1920,2), end=c(1939,2)) 
CornF.ts = ts(data=hamilton_final_CornF, frequency = 3, start=c(1920,2), end=c(1939,2)) 
OatsS.ts = ts(data=hamilton_final_OatsS, frequency = 3, start=c(1920,2), end=c(1939,2)) 
OatsF.ts = ts(data=hamilton_final_OatsF, frequency = 3, start=c(1920,2), end=c(1939,2)) 
RyeS.ts = ts(data=hamilton_final_RyeS, frequency = 3, start=c(1920,2), end=c(1939,2)) 
RyeF.ts = ts(data=hamilton_final_RyeF, frequency = 3, start=c(1920,2), end=c(1939,2)) 
LardS.ts = ts(data=hamilton_final_LardS, frequency = 3, start=c(1920,2), end=c(1939,2)) 
LardF.ts = ts(data=hamilton_final_LardF, frequency = 3, start=c(1920,2), end=c(1939,2)) 
TBill.ts = ts(data=hamilton_final_TBill, frequency = 3, start=c(1920,2), end=c(1939,2)) 
# Calculating Expected Commodity Price Inflation an time = t
CornPred.ts <- 100 * (log (CornF.ts) - log(CornS.ts))
OatsPred.ts <- 100 * (log(OatsF.ts) - log(OatsS.ts))
RyePred.ts <- 100 * (log(RyeF.ts) - log(RyeS.ts))
LardPred.ts <- 100 * (log(LardF.ts) -log(LardS.ts))

# Looking at seasonality of the Expected differences at t
#Corn first
SeasC1 <- cycle(CornPred.ts)
CornPred.lm <- lm(CornPred.ts ~ factor(SeasC1))
summary(CornPred.lm)
# names(CornPred.lm)

coef(CornPred.lm)
coef(summary(CornPred.lm))[,2] # column 2 is standard error. Could have also written [, "Std. Error"]
anova(CornPred.lm)

# Now the others
SeasC2 <- cycle(OatsPred.ts)
OatsPred.lm <- lm(OatsPred.ts ~ factor(SeasC2))
summary(OatsPred.lm)

# These rye numbers are slightly off
SeasC3 <- cycle(RyePred.ts)
RyePred.lm <- lm(RyePred.ts ~ factor(SeasC3))
summary(RyePred.lm)

SeasC4 <- cycle(LardPred.ts)
LardPred.lm <- lm(LardPred.ts ~ factor(SeasC4))
summary(LardPred.lm)

# Looking at the seasonality of the actual differences at t
LagCornS.ts <- lag(CornS.ts, k=1)
LagOatsS.ts <- lag(OatsS.ts, k=1)
LagRyeS.ts <- lag(RyeS.ts, k=1)
LagLardS.ts <- lag(LardS.ts, k=1)
CornSSDiff.ts <- 100 * (log (LagCornS.ts) - log (CornS.ts))
OatsSSDiff.ts <- 100 * (log (LagOatsS.ts) - log (OatsS.ts))
RyeSSDiff.ts <- 100 * (log (LagRyeS.ts) - log (RyeS.ts))   
LardSSDiff.ts <- 100 * (log (LagLardS.ts) - log (LardS.ts))  
                             
SeasCC1 <- cycle(CornSSDiff.ts)
CornSSDiff.lm <- lm(CornSSDiff.ts ~ factor(SeasCC1))
summary(CornSSDiff.lm)

SeasCC2 <- cycle(OatsSSDiff.ts)
OatsSSDiff.lm <- lm(OatsSSDiff.ts ~ factor(SeasCC2))
summary(OatsSSDiff.lm)

SeasCC3 <- cycle(RyeSSDiff.ts)
RyeSSDiff.lm <- lm(RyeSSDiff.ts ~ factor(SeasCC3))
summary(RyeSSDiff.lm)

SeasCC4 <- cycle(LardSSDiff.ts)
LardSSDiff.lm <- lm(LardSSDiff.ts ~ factor(SeasCC4))
summary(LardSSDiff.lm)

# To extract coefficients "Estimate" and SEs "Std. Error"
reportedcoeffN <- function(modelobject,nn,what){
  g <- summary(modelobject)$coefficients[[nn,what]]
  return(g)
}

# To extract F p-values 
lmp <- function (modelobject) {
  if (class(modelobject) != "lm") stop("Not an object of class 'lm' ")
  f <- summary(modelobject)$fstatistic
  p <- pf(f[1],f[2],f[3],lower.tail=F)
  attributes(p) <- NULL
  return(p)
}
# To extract F values
fvlu <- function (modelobject) {
  f <- summary(modelobject)$fstatistic[["value"]]
  return(f)
}
# Testing the above functions

reportedcoeffN(LardPred.lm,2,"Estimate")
fvlu(LardPred.lm)
lmp(LardPred.lm)



sink("~/Dropbox/Replication Workshop/Student Papers/Saleuddin/Analysis/R/table3.out")
# Replicate Table 3
"Table 3 - Regressions of Expected and Actual Commodity Price Inflation on Constant and Seasonal Dummies"
"Corn"
Table3C <- matrix (c(reportedcoeffN(CornPred.lm,1,"Estimate"),reportedcoeffN(CornPred.lm,3,"Estimate"),reportedcoeffN(CornPred.lm,2,"Estimate"),fvlu(CornPred.lm),reportedcoeffN(CornPred.lm,1,"Std. Error"),reportedcoeffN(CornPred.lm,3,"Std. Error"),reportedcoeffN(CornPred.lm,2,"Std. Error"), lmp(CornPred.lm),reportedcoeffN(CornSSDiff.lm,1,"Estimate"),reportedcoeffN(CornSSDiff.lm,3,"Estimate"),reportedcoeffN(CornSSDiff.lm,2,"Estimate"),fvlu(CornSSDiff.lm),reportedcoeffN(CornSSDiff.lm,1,"Std. Error"),reportedcoeffN(CornSSDiff.lm,3,"Std. Error"),reportedcoeffN(CornSSDiff.lm,2,"Std. Error"), lmp(CornSSDiff.lm)), ncol=4, byrow=TRUE)
colnames(Table3C) <- c("Constant", "September dummy", "May dummy", "F test and F test P value")
rownames(Table3C) <-c("Expected Inflation: Coeff", "Expected Inflation:SE", "Actual Inflation: Coeff","Actual Inflation: SE")
Table3C <- as.table(Table3C)
Table3C


"Oats"
Table3O <- matrix (c(reportedcoeffN(OatsPred.lm,1,"Estimate"),reportedcoeffN(OatsPred.lm,3,"Estimate"),reportedcoeffN(OatsPred.lm,2,"Estimate"),fvlu(OatsPred.lm),reportedcoeffN(OatsPred.lm,1,"Std. Error"),reportedcoeffN(OatsPred.lm,3,"Std. Error"),reportedcoeffN(OatsPred.lm,2,"Std. Error"), lmp(OatsPred.lm),reportedcoeffN(OatsSSDiff.lm,1,"Estimate"),reportedcoeffN(OatsSSDiff.lm,3,"Estimate"),reportedcoeffN(OatsSSDiff.lm,2,"Estimate"),fvlu(OatsSSDiff.lm),reportedcoeffN(OatsSSDiff.lm,1,"Std. Error"),reportedcoeffN(OatsSSDiff.lm,3,"Std. Error"),reportedcoeffN(OatsSSDiff.lm,2,"Std. Error"), lmp(OatsSSDiff.lm)), ncol=4, byrow=TRUE)
colnames(Table3O) <- c("Constant", "September dummy", "May dummy", "F test and F test P value")
rownames(Table3O) <-c("Expected Inflation: Coeff", "Expected Inflation:SE", "Actual Inflation: Coeff","Actual Inflation: SE")
Table3O <- as.table(Table3O)
Table3O

"Rye"
Table3R <- matrix (c(reportedcoeffN(RyePred.lm,1,"Estimate"),reportedcoeffN(RyePred.lm,3,"Estimate"),reportedcoeffN(RyePred.lm,2,"Estimate"),fvlu(RyePred.lm),reportedcoeffN(RyePred.lm,1,"Std. Error"),reportedcoeffN(RyePred.lm,3,"Std. Error"),reportedcoeffN(RyePred.lm,2,"Std. Error"), lmp(RyePred.lm),reportedcoeffN(RyeSSDiff.lm,1,"Estimate"),reportedcoeffN(RyeSSDiff.lm,3,"Estimate"),reportedcoeffN(RyeSSDiff.lm,2,"Estimate"),fvlu(RyeSSDiff.lm),reportedcoeffN(RyeSSDiff.lm,1,"Std. Error"),reportedcoeffN(RyeSSDiff.lm,3,"Std. Error"),reportedcoeffN(RyeSSDiff.lm,2,"Std. Error"), lmp(RyeSSDiff.lm)), ncol=4, byrow=TRUE)
colnames(Table3R) <- c("Constant", "September dummy", "May dummy", "F test and F test P value")
rownames(Table3R) <-c("Expected Inflation: Coeff", "Expected Inflation:SE", "Actual Inflation: Coeff","Actual Inflation: SE")
Table3R <- as.table(Table3R)
Table3R

"Lard"
Table3L <- matrix (c(reportedcoeffN(LardPred.lm,1,"Estimate"),reportedcoeffN(LardPred.lm,3,"Estimate"),reportedcoeffN(LardPred.lm,2,"Estimate"),fvlu(LardPred.lm),reportedcoeffN(LardPred.lm,1,"Std. Error"),reportedcoeffN(LardPred.lm,3,"Std. Error"),reportedcoeffN(LardPred.lm,2,"Std. Error"), lmp(LardPred.lm),reportedcoeffN(LardSSDiff.lm,1,"Estimate"),reportedcoeffN(LardSSDiff.lm,3,"Estimate"),reportedcoeffN(LardSSDiff.lm,2,"Estimate"),fvlu(LardSSDiff.lm),reportedcoeffN(LardSSDiff.lm,1,"Std. Error"),reportedcoeffN(LardSSDiff.lm,3,"Std. Error"),reportedcoeffN(LardSSDiff.lm,2,"Std. Error"), lmp(LardSSDiff.lm)), ncol=4, byrow=TRUE)
colnames(Table3L) <- c("Constant", "September dummy", "May dummy", "F test and F test P value")
rownames(Table3L) <-c("Expected Inflation: Coeff", "Expected Inflation:SE", "Actual Inflation: Coeff","Actual Inflation: SE")
Table3L <- as.table(Table3L)
Table3L




sink()


